Equivalence of operator-splitting schemes for the integration of the Langevin equation 
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We investigate the equivalence of different operator-splitting scfiemes for tfie integration of the 
Langevin equation. We consider a specific problem, so called the directed percolation process, which 
can be extended to a wider class of problems. We first give a compact mathematical description 
of the operator-splitting method and introduce two typical splitting schemes that will be useful in 
numerical studies. We show that the two schemes are essentially equivalent through the map that 
turns out to be an automorphism. An associated equivalent class of operator-splitting integrations 
is also defined by generalizing the specified equivalence. 

PACS numbers: 



Several kinds of models for lattice-based dynamic pro- 
cesses have been studied for decades to understand the 
characteristics of non-equilibrium systems, especially, fo- 
cused on the critical phenomena. To mention a few, di- 
rected percolation (DP) Q, S H |1 H 111 , contact pro- 
cess 't',^, and catalytic reactions^ are such examples. 
In these studies, it was found that, in spite of diversity in 
microscopic details, various models exhibit critical phe- 
nomena that are essentially identical to those for the ab- 
sorbing phase transition in the DP model and the DP 
universality class was inferred for them 0, ^| . 

In order to explain the DP universality class from the 
unified point of view, so called the DP Langevin equation 
is proposed 0,0] as 



dtp = ap- bp^ + D\/^p + a^p 



(1) 



where p = p(r, t) is a non-negative field variable for con- 
centration and ^ is a white noise with zero mean satisfy- 
ing (^(r,i)C(r',t')) = 2(5(r - r')S(t - t'). The coefficient 
a is the tuning parameter for the phase transition, and 
6, D, a are positive constants. 

The field theoretical approach to Eq. named as 
the Reggeon field theory 12], unveils the critical behavior 
of the DP universality class for spatial dimensions com- 
parable to or higher than the critical one dc = 4. The 
mean field theory applies in high dimensions {d > dc ~ 4) 
and the perturbative results are well established near and 
lower than dc through the standard e = dc — d expan- 
sion. In lower dimensions, the series expansion meth- 
ods for lattice-based DP models provide the most precise 
estimates for the critical scaling exponents 

111, 

which 

have been also confirmed by extensive numerical simula- 
tions [loj . 

Attention has also been paid to the direct numerical in- 
tegration of Eq. , especially, in the quantitative study 
of the lower-dimensional cases This seems to be 

a simple numerical integration of the partial differential 
equation at a glance. However, as long as the absorb- 
ing transition is concerned, one has to meet an annoying 
block which is by no means easily tractable. The conven- 



tional Euler integration technique using a discrete time 
interval At may result in a negative value of p due to 
the uncontrolled random noise term, and then any fur- 
ther sensible integration of the equation is impossible. 
In particular, this nuisance may appear very easily for 
small p where the noise term (~ ^/p) is dominant over 
the deterministic term (~ p). As the absorbing critical 
behavior occurs in the p — ^ limit, the proper treatment 
to guarantee a small but positive value of p is not only 
of a technical interest, but also a critical issue in the nu- 
merical study of the absorbing phase transition. Some 
numerical schemes have been put forward to overcome 
such numerical fragility d4|J , but without much success. 

Recently, Dornic et al. 0| utilized the operator- 
splitting method for integrating the Langevin equations 
(similar idea earlier in Ref. describing various kinds 
of absorbing critical phenomena 0| . In this method, the 
evolution in time is divided into two parts, each of which 
is exactly solvable. The successive integration of the two 
parts during At is regarded as one-step time evolution 
of Eq. (QJ, getting exact as At 0. This method has 
a couple of outstanding advantages over the preceding 
ones [l^. As the noise term is treated exactly, the non- 
negativity of p is always guaranteed. Therefore one can 
use a relatively large value for At to stay in the critical 
region where p is small, which can also save the com- 
puting time considerably. An astonishing observation is 
that the critical behavior seems to be fairly insensitive to 
the magnitude of At. In fact, even if At = 0.25 in their 
paper, the critical point is shifted only by 1% from the 
extrapolated value in the limit of At — > 0. They have 
reported many successful applications of this method in 
various types of absorbing critical phenomena. 

It is a tough task to rigorously explain why the 
operator-splitting method with a relatively large time 
interval yields a reliable result for the critical behavior 
of a certain problem and find what is the criterion for 
this approach to be valid. Our preliminary work based 
on the perturbative expansion in At reveals that the 
operator-splitting scheme renormalizes the given coeffi- 
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cients (a, b, D, a), generates higher order terms, and 
new non-Gaussian noises. Furthermore, we find that it 
may be possible for the renormahzed and newly gener- 
ated coefficients to change their signs, which could pro- 
duce a diverging solution, a first-order phase transition, 
or a higher-order multicritical point. Unfortunately the 
coefficients are expressed in terms of alternating infinite 
series, from which it seems hardly probable to derive the 
validity criteria of At for maintaining the DP critical be- 
havior. 

In this paper, we present a compact mathematical 
description of the operator-splitting method. We ana- 
lytically show that some of seemingly different splitting 
schemes are mathematically equivalent in the sense that 
there is an exact transformation map relating different 
splitting schemes. We hope that our result may eluci- 
date the structure of the At-dependent terms and even- 
tually help us understand the characteristics of the upper 
bound for At in general operator-splitting schemes. 

First, we summarize the operator-splitting scheme 
specified in Ref. and then provide a compact math- 
ematical description for it. We next consider another 
scheme with a different choice in splitting the dy- 
namic process. Using the Baker-Campbell-Hausdorff 
(BCH) formula, we prove the equivalence of two differ- 
ent schemes with the analytic expression of the parameter 
transformation function. 

In the numerical study, the embedding space is re- 
placed with a mesh-like lattice with lattice constant Ax. 
Then, at site i, Eq. becomes an ordinary differential 
equation for pi{t) 

= ap, - bpf + ^ e^jPj + (2) 

3 

where pi is the time derivative of pi, a = a — 2dD/ (Ax)^, 
eij = 1 for nearest neighbor sites i and j or eij = 
otherwise, and d is the spatial dimension. Eq. is 
a set of coupled equations where the dynamics at site 
i is influenced by field variables pj at nearest neigh- 
bors. As an additional approximation, we assume that 
Pj is piecewise constant with an initial value pj (t) during 
the integration from t to t + At. During this interval, 
Eq. ^ then becomes decoupled with an effective field 
Ci = D/ [AxY X]j ^ijPjit)^ which leads to 

p = c + dp- bp"^ +(T^ (3) 

where the subscript i is dropped for simplicity. 

The main idea of the operator-splitting method is to 
separate the right hand side of Eq. Q into two parts, 
each of which can be treated exactly. In particular, it 
is important to have or find the exact solution of the 
Fokker-Planck (FP) equation associated with the part in- 
cluding the noise term, guaranteeing the non-negativity 
of p. Dornic et al. [31 considered a stochastic equation 



without the nonlinear term and a purely determinis- 
tic equation involving the nonlinear term only, given as 

p^c + dp + i , p^ -bp^ . (4) 

The FP equation associated with the first stochastic 
equation can be solved exactly 0, . The conditional 
probability density Vsip,t] po) can be obtained analyti- 
cally for an initial value po, vanishing for p < at any 
time. The deterministic equation is trivially integrated 
as pD(t; Po), which also preserves the non- negativity of p. 

The numerical integration during At is done as follows: 
Given an initial value pt at time t, p is updated first 
by sampling a value appropriate to the exact probability 
distribution Vs{p, At; pt), which is denoted by ps(Ai; pt). 
Then, resetting this updated value as an initial value at 
time t, we integrate the deterministic equation over At, 
which yields 

Pt+At = pu{At; ps{At; Pt)) . (5) 

The same procedure is performed for all sites in parallel. 
This constitutes a single step in discrete time dynamics 
with the time interval At. The next step follows with 
new initial values of pt+At and newly tuned values of c. 
This is the idea introduced in Ref. |l5| . 

The splitting scheme can be described mathemati- 
cally in terms of the probability density 7^(p, t). The 
probability density after a single step can be written 
by P(p,t + At) = C{a,b,D,a)V{p,t), where C is the 
time evolution operator. Notice that C is the product 
of two consecutive evolution operators, having the form 
of e^*^De^*^s^ where Ls.d are the Fokker-Planck (FP) 
operators |0| associated with the stochastic and deter- 
ministic differential equations, respectively in Eq. |0J. In 
the Ito calculus, one writes 

£ ^ gAtiogAtLs ^ gAiPf,p2gAt(-P(c+ap) + PV=^p) ^ ^g-j 

where P = d/dp. We remark that Eq. ijHJl compactly 
contains the whole information of the operator-splitting 
method represented by Eqs. and 0. 

The exact FP operator is given by L = Ls + -^d • Due 
to the non-cummutativity of two operators, [Lq, Ld] 7^ 0, 
the exact evolution operator >Coxact = e"^*^ differs from 
the operator-splitting evolution operator C in higher or- 
ders of At. The difference between >Coxact and C can be 
found systematically in power series of At using the BCH 
formula: for = e^e^ , Z can be expressed as Ho] 

Z = A+ [ dtg {e^'^^ e* '"^'') B , (7) 
Jo 

where g{x) = 1 + Y.Z=i ''r^{l+i) " 1)"" and adx is a 
linear map of which the operation is defined by adxY = 
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[X, Y] . This leads to the rather famihar BCH formula; 

Z = A + B+^[A,B] + j^[A,[A,B]]-j^[B,[A,B]] 
+ .. + k^,,..,^Jwi,[..[wn,[A,B]]..]] + ... , (8) 

where Wi stands for either A or _B, and A:m)i,..,u)„ is a 
constant scalar. Here we do not write down the explicit 
expression of kn]i,..,w„, but only note that its sign con- 
stantly changes with n. 

With A = AtLo and B = AtLs as in Eq. ©, the 
commutator [A, B] produces a new type of higher-order 
noise terms such as P^p^ by the interplay of the non- 
linear term Pp^ and the noise term P^p. Through the 
nested commutators in Eq. ((SJ, the operator-splitting FP 
operator InC = Z includes not only higher order deter- 
ministic terms like Pp" but also higher order noise terms 
like P"^p". The coefficients of the preexisting lower order 
terms such as Pp, Pp^, and P^p are also modified. In 
the renormalization group sense, the higher order terms 
are usually irrelevant to the DP critical behavior, but 
only when the appropriate stability condition is satisfied. 
For example, the fixed point solution of the determinis- 
tic part in the operator-splitting FP operator should not 
have a different structure from that in the exact FP op- 
erator. Therefore the stability condition depends criti- 
cally on the detailed Ai-dependence of the coefficients. 
Unfortunately, the complexity of fc^j, do not allow 
us to judge the stability criteria in a sensible way. Any 
conclusion derived from a truncated finite series in the 
perturbative expansion of Eq. Q may not contain rel- 
evant information on the stability, especially due to the 
alternating nature of the series. It seems also impossible 
to sum up the infinite series in a closed form even for 
the coefficients of the lower order terms. Thus it is not 
clear that the modified coefficients due to the operator- 
splitting method employed in Ref. jl^ still guarantee the 
stability of the DP-type solutions for any value of At. 

Now we focus on classifying various operator-splitting 
schemes into equivalent classes, which will greatly re- 
duce the efforts to derive the stability criteria for At 
in general operator-splitting methods. First, we notice 
that the nested commutators in Eqs. (TJ and © can 
be easily summed up when the commutator [A, B] can 
be written as a linear combination of A and B only. 
Consider the simple case of [A, B] — -fB with a real 
constant 7. As ad^i? 0, it is easy to show that 

^(gad.gtad«)5 - (1+E™=1 StTIT^^^-I)")^- ^^^^^ 

we obtain 

Z^A + a^B, (9) 

where = 7/(1 — e^^). Therefore we find e^e^ — 
^A+a^B equivalently e^e^ — e<^-t^+^_ A general case 
of [A, B] = ^aA + ^bB can be reduced to the simple case 
by replacing B hy B proportional to the commutator. 
The rather complicated result is not shown here. 



Next, we observe that the components of the FP op- 
erator L satisfy the above special commutation relation 
such as [Pp, P] = -P, [Pp, Pp^] = Pp^, and [Pp, P'^p] = 
—P^p. This implies that the linear term Pp can move 
around rather freely between Ls and Ld without caus- 
ing too much complication in the operator-splitting FP 
operator Z . We note that the other commutators do not 
satisfy the special relation and the algebra is not closed. 

Consider a different splitting scheme where the linear 
term is included in the deterministic part. The splitted 
Langevin equations are 

/5 = c + cr^^, p^ap-bp^ . (10) 

As in the previous splitting in Eq. I^J , both equations can 
be treated exactly. The corresponding evolution operator 
can be written as 

C'{a,b,D,a) = e.''tP(bp-^--p)e.^ti-Pc+P'-'p), (n) 

Using the identity in Eq. © , we split the first exponential 
map as 

L'{a,b,D,a) = e^*^(^/"^^*)'''e-^*^"''e^*(-^^+^'"''') . 

(12) 

The last two exponential maps can be merged together 
as 

(13) 

By comparing Eqs. JBJ and (|13() . one establishes the 
relation between the two operator-splitting schemes as 

C'{a, b, D, a) = C{a' , b' , D' , a') , (14) 

where 

a' = a + 2dD{aaAt - l)/(Aa;)2 , b' = 6/aaAt , 
D' = aaAt.D , a' = V^aAtC • (15) 

Since a is positive-definite and does not vanish nor di- 
verge for any finite aAt, the transformation between C 
and £' forms an automorphism (one-to-one and onto it- 
self) in the parameter space of {a,b, D, a). The trans- 
formation preserves the sign of the parameters (5, D, a) 
except a (tuning parameter) where the reformulation of 
the discrete Laplacian is involved. Therefore, the differ- 
ent operator-splitting methods are related to each other 
only by trivial rescaling of parameters with a shift of the 
critical point. Any phenomenon observed in one split- 
ting scheme is also expected in the other scheme and 
the stability conditions for At can be traced using the 
transformation of Eq. (|15(l if it is known for one specific 
operator-splitting method. 

We may consider a more general splitting where the 
linear term is arbitrarily divided into two parts. That is, 
the splitted Langevin equations are now 

p — c+ {a ~ k)p + Gyfp ^ , p — kp—bp^. (16) 
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where fc is a real arbitrary constant. Note that we can 
still integrate both equations exactly. The corresponding 
time evolution operator is given as 

Ckia,b,D,a) = e'^tPibp--kp)^At{-Pic+ia-k)p)+p-a^p)_ 

(17) 

Similarly, we find 

(18) 

where 

- goAt _ gfeAt 

Hence one can generalize Eqs. H14f) and (|15|l as follows; 

Ck{a,b,D,a) = £{ak,bk,Dk,ak) , (20) 

where 

flfc = a + 2rfi:'(/3fe - l)/(Aa:;)2 , bk = b/okAt , 
Dk = PkD , au = . (21) 

Due to the property of (3, inherited from that of a, the 
transformation between Ck and C also forms an automor- 
phism for any k. Consequently, the solution-structure 
yielded by Ck is always preserved irrespective of k, and 
thus such operations can be represented by C This di- 
rectly demonstrates that Ct 's form an equivalent class of 
operator-splitting integration of the DP Langevin equa- 
tion. 

In summary, we present a compact mathematical 
description of the so-called operator-splitting method, 
which was proposed in Refs. |l5[ ll(il | for the numeri- 
cal integration of the DP Langevin equation of Eq. 
Based on this, we show analytically that some splitting 
methods are mathematically equivalent with the explicit 
transformation function of the model parameters. Con- 
sequently, we find that the splitting methods Ck^ form 
an equivalent class of integration in the sense that the 
solution-structure and the transition property between 
the solutions are always conserved for any k. In the 
meantime, we also address that the difference between 
the original dynamics and that by the operator-splitting 



scheme is traceable by the perturbation theory of the pre- 
sented mathematical description. However, the informa- 
tion from the perturbation theory seems not sufficient to 
decide whether the splitting scheme still preserves the es- 
sential feature of the original dynamics. Nonetheless, our 
work on the equivalence class will be of considerable help 
to examine the validity of the operator-splitting scheme 
in studying the universality of the DP Langivin equation. 
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